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io INTRODUCTION 


ee Se ee ee ee ee eee 


A. BACKGROUND 


A Magnetchydrodynamic (MHD) generator is a device for 
Menverting the energy of a flowing gas, or of a liguid 
metal, directly into electrical energy. As With 
conventional electrical generators, the conducting medium 
crosses a magnetic field; however, MHD deals with conducting 
piuids instead of solid moving parts. Anticipated 
efficiencies are far above those of conventional generators 
because of the higher temperatures in the conversion 
Channel. The thermal energy of flowing gaseS is changed 
directly to electrical energy, eliminating the mechanical 
energy step cf the conventional generator which means that 


there are fewer parts to wear out. 


Although MHD power generation has been studied for 
Nearly 30 years, financing of MHD research waned in the 
1960's. Interest is now being rekindled because of a new 
"energy consciousness" of government, industry, and the 
fmeeeitary. MHD offers an attractive alternative energy 
transformation means which is relevant to the pursuit of 
clean ways of using the vast coal reserves in the United 
States and waking more efficient use of dwindling petroleum 
resources. MHD promises to fulfill the requirements for a 
light-weight, high-power source for military use. Hor 
eXample, the Navy, recognizing itsS energy requirements and 
dependency on world fossil-fuel supplies{1], 1S sponsoring 
research on practicai MHD devices with large power outputs 
and high efficiencies[2]. The Air Force is interested in 


power sources for airborn weapons systems. 


The basic MHD device consists of a channel through which 
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ionized jas or liquid netal flows. A magnetic field is 
imposed perpendicular to the flow direction; as a result, an 
electric field is induced and power 1s tapped by means of 
electrodes. Principal concern 1s with ionized gas 
(hereafter refered to as plasma) deviceS since with these 
the combustion energy can be more directly changed to 


electrical power. 


A major problem involved in the performance of MHD 
generators 1s the high voltage loss in the vicinity of the 
electrodes. Accurate performance predictions in tne design 
of MHD devices require a realistic determination of these 
voltage losses. This research investigates the nature of 
electrode voltage drops and computer models are developed to 
describe the physical phenomena involved. With the results 
from this investigation, appropriate steps may be taken to 


Mentmize the losses. 


The principal voltage loss mechanisms in the MAD 
generator can be divided into two main classes, ohmic and 
Sheath losses. Ohmic drops are those that occur because of 
the finite conductivity of a real plasma. Thermal boundar 
layers, degree and kinetics of lonization, and Joule heating 
are factors affecting the ohmic resistivity of the plasma. 
Sheath drcps occur as the result of the Debye shielding 
which forms a non-neutral layer adjacent to the electrode 
and results ina space charge field. Material prcblems 
restrict the temperatures at which the electrodes can 
Operate. In many cases cooling of the electrodes is 
required since the plasma, in order to maintain a high 
lonization, must be hotter than the working temperatures of 
most materials. This temperature difference between the 
electrodes and the plasma further aggravates the voltage 
losses because of the presence of the thermal boundary 
layers. As will be shown, voltages losses can be as much as 


50% or more of the total power output, with sheath drops 
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accounting for a non-negligible fraction of the drop and 


boundary layer losses the rest. 
PameeneVIEN OF PREVIOUS WORK 


Analyses and experiments on boundary layer phenomena in 
MHD generators are numerous. It 1s not the purpose of this 
work to further add to the literature on the boundary layer 
messes, but, aS wiil be explained in Part C of this section, 
this phase of the problem was treated aS a necessary 
appendage to the problem of the sheath drop. Consequently a 
short review cf boundary layer work is presented here: High 
and Felderman[ 3] as well as Argyropoulos, Demetriades, and 
Lackner{4]} studied the problem of the turbulent boundary 
layer and several descriptions of the nature of the 
temperature field. Doss, Dwyer and Hoffman[5} investigated 
the boundary layer ana used a rudimentary collisionless 
Sieath as an "inside" boundary condition. Kessler and 
Eustis[6] rerorted on experiments with electrode temperature 
effects on turbulent boundary layers, and Rubin and 
Eustis[ 7} extended Kessler's work to include the effects of 
electrode size. Wu et al.{38} and Oliver and Mitchner[ 9] 
investigated non-uniformities of current Gres iede ot OS 


Within the boundary layer. 


Although the existence of the sheath is well understood, 
its effects have been investigated toa a much lesser extent 
than these of the boundary Jlayer. fies 1S  pamime ip 4) iy 
because this loss is described by a relatively complicated 
Set Of cCoupied, non-linear, partial differential equations 
that present considerable deeteere ult y for numerical 
solutions. There is no known exact solution for these 


eguations. 


Most work on sheath phenomena is embodied in "pgroke" 


theory i: estigations, that is in the mutual effects on a 
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quiescent plasma which is disturbed by the presence of an 
electric prcbe. Such work is relevant to MHD electrodes 
Since the anode is essentially a heavily biased probe in 
contact with a plasma, which is quiescent within the sheath 
region (see Part C of this section). The cathcde, which is 
an electron emitter, has some of the Same characteristics as 
the anode, but is more complicated £0 represent 
analytically. The description of the cathode will not be 


attempted here. 


Necessary simplifications limit existing probe solutions 
to special cases as can be seen by sampling tne literature 
regarding plasma probes. Lam{10] solves a sheath problem hy 
matching the boundaries of "“anner" and "outer" solutions and 
discusses their dimensionality. Stahi and Su{ 11] use the 
Same approach of separating the sheath, ambipolar region and 
free strean. Additionally, they prove the existence of a 
emeath on a flat probe. Cohen{ 12] criticizes this approach 
Or Separate regions on mathematical grounds saying that 
gGuantities are forced to fit and some derivatives tend to be 
discontinuous. lickKee and Mitchnerf 13} deal with a 


collisionless sheath (zr ) Duteerncluce 2Onization “ana 
S 


recombination 16m the ambipolar reqlon. Balley and 
Touryan{ 14] investigate a sheath that is large enowuh to be 
of the same order of magnitude as the boundary layer, and 
take advantage of Blasius' similitude co-ordinates to reduce 


the problem to one dimension. 


Several solutions are available for spherical probes 
ancluding Kiel{15] , Barad and Cohen[ 16} , Su and Lanm[17]} , 
and Cohen{12]. Kiel ignored diffusion while Barad and Cohen 
neglected ionization and recombination. Su and Lam, ana 
Cohen were the first to use a systematic analysis of probes 
in collisicn-dominated plasmas. Lengyel[{18], by dropping 


diffusion, was able to modify the elliptic equations in such 
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a way as to solve them by the method of characteristics as 


one might do with hyperbolic equations. 


Gung, Talbot, and Touryan[ 19}, in a rather extensive 
review of probe work, state that no general soluticn is 
available EOE determining charge density and svecies 
temperature for probes smali relative to boundary layer 
thicknesses. The work of this thesis may help to fill that 
void because the model of the anode is that of an array of 
point or line probes immersed in the non-convective portion 
of a boundary iayer. Results of this work nay also be 


useful to other fields such as arc-discharges and lasers. 


Meee PHYSICAL MODELING 


Both theoretical and experimental results indicate that 
in the MHD environment, where pressures are near 
atmospheric, the Sheath lies well within the boundary layer 
(about 10-5 m for the sheath thickness or three orders of 
magnitude less than the boundary layer). Additionally, the 
Sheath is only about one or two orders of magnitude larger 
than the electron mean-free-path for combustion gases. It 
becomes apparent then, that convection plays little cr ne 
role in the behavior of the Sheath other than to modify the 


gross temperatures encountered near tne wall. Conseguentiy, 
the boundary layer or ohmic problem can be divorced from the 
Sheath and treated separately[20]. The boundary layer, as 
well as the wall conditions, determines the gaS tenperature 
for the Sheath region. Thus, the boundary layer voltage 
drop can be added to the sheath drop to give the total loss 


due to electrode effects. 
As mentioned earlier, analyses of boundary layer effects 


in MHD have keen developed, but they tend tc be conplicated 


and difficult to use. In order to have a complete picture 
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of electrode drops and to assist in interpreting data, a 
Simplified method for determining the voltage drop across a 
thermal boundary layer haS been devised. It has the 
Capability cf rendering quick results with little computer 


Space and time. It does not pretend to be a _ complete 


analysis, but through appropriate assumptions and 
Simplifications it gives the drop as a fe yaver calle yy of 
temperature only. Appendix A gives the details of the 


method along with its limitations. 


An investigation of the nature of the controlling 
eguations for the sheath shows that a one-dimensional 
solution does not exist except for very special cases. 
This investigation resulted in the publication of Ref. 21, 
and a summary of that paper may be found in Appendix EE. {It 
provides a basis for many of the assumptions used in this 
work. The following concluSions which appear in Ref. 21 
helped to shape the model chosen for this9 wock: Cucrent 
density must decrease away from the electrode by whatever 
means possible. The mechanisms of geometry, current 


constrictions at the electrode, or ionizaticn/recombination 


can effect this decrease Singly OE FeCl Clay. A 
two-dimensional flat plate in the absence or 
ijonization/recombination would be at least partially 


One-dimensional away from the edges and is therefore not a 


workable model. 


Since previous solutions have been obtained with 
Spherical probes, it appears that a geometrical decrease of 
current density 1S a proper means of satisfying the nature 
of the equations. Chen{22] used ionization/recombination to 
extract a strictly one-dimensional solution for the sneath. 
What remains then is to study the features that result fron 
current constrictions at discrete periodic sites rather than 


uniforn current density across the electrode surface. 
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PeeeeOoJECTIVES OF THE PRESENT WORK 


In this WOLK, a two-dimensional, Elvat plate, 
non-emitting electrode .is modeled. Here, the current has 
only discrete points through which to flow on the electrode 
Surface (See Fig. 1). These points are assumed to be 
periodic. The points simulate the current constriction at 
the electrode which may be there because of temperature 
irregularities or surface roughness. A collection of these 
points or nodes might represent the anode spots observed by 
Beeolarz{ 23} and Kimblin[ 24 }. In a two-dimensional 
representaticn the point is actually a "wire" of unit depth. 
The Spacing cf the active sites must be greater than the 
sheath thickness in order to avoid ore-dimensional effects, 
and small enough to be compatible with a two-dinensional 
computational field. It must be emphasized here that this 
work does not consider the phenomenon of arc discharge, 
which is the result of thermal instabilites[{25]. kather it 
Mens at the current constriction required to satisfy 


memtinuity of charge, Ohm's law, and Poisson's equation. 


The primary objective of this work is to present a 
mathematical model of a two-dimensional (Cartesian ceometry) 
electrode in a quiescent MHD plasma. The COnErolLling 
equations are solved for the potential, the current, and the 
Charge density distributions in the vicinity of the 
electrode. The additional effects of a magnetic field and 
Joule heating are investigated. It appears that this is the 
first investigation of the magnetic field within the sheath. 
The resulting program generates sheath and ambipolar regicns 
in a self-consistent way using the same set of equations 
throughout the field, obviating the need to match layers. 
The size ot the sheath and the voltage drop attributable to 


its shielding effects are investigated. 
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E. ORGANIZATION 


In Section II the physical concepts and the controiling 
and boundary equations are developed. The equations are 
cast into a form for use in the computer in Section III. [ft 
includes a description of the numerical methods which 
meogauveed Successful results, as well as a short chronology 
of techniques that failed to work. A summary of the results 
of potential and charge distributions, and current and 
temperature plots for all the cases studied 1S included in 
Section IV. Section V lists the conclusions reached fron 
both a physical and numerical point of view and gives 


examples on the use of the results. 


Appendix A shows the technigue for determining the 
boundary layer losses. Appendix B is a condensaticn of 
Ref. 21 proving the non-existence of a one-dimensional 
solution for the limit of the non-reacting flux of charges 
in a collision-dominated plasma. A dimenSional analysis 
which derives a set of independent non-dimenSional variables 
moe the controlling equations is found in Appendix Cc. 
Metaohnally, Appendix C contains a fractional analysis of 
these and other equations to estimate the significance of 
certain physical effects. Finally, the computer prograns 
used in this work are included in the section "Conputer 


Programs", with explanations as to their uSe. 
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A. THE SHEATH ENVIRONMENT 


A sheath is a non-neutral region which lies adjacent to 
the electrode or insulator surface. The strict definition 
of a plasma requires charge neutrality, thus the sheath is 
mete part of the plasma since a space charge exists. For 
example, the anode is positive, which means that it attracts 
electrons and rerels positive ions. As the electrons 
collect at the anode surface, the anode potential is 
partiaily shielded from the rest of the plasma by the 


space-charge potential drop. 


Between the sheath and the free stream lies the 
aMbipolar region. As the name implies, it has an eauval 
number of positive and negative charges. [In this transition 
region the electrons are slowed by heavier ions, ard strona 
concentration gradients exist. The electrodes themselves 
are usually metallic, and may be pins, wires, or flat plates 
separated by ceramic insulator wall segments (segmented 
electrodes). The surface or the electrodes, even when highly 
polished have roughnesses of the order of the sheath 
thickness, and in addition exchange material with the 
meowing plasma, further adding to the irregularities. These 
irregularities increase the tendency for the current to 
constrict by providing active sites for the current to flow 


along minimum energy paths. 


The size of the sheath in #WHD generators is about 
10-5 m, putting it well within the boundary layer whicn is 
about 10-2 mm. The flow is therefore essentiaily stagnated 
in the sheath regicn. Though the ambipoiar region extends 


further into the boundary layer, a fractional analysis 
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(Appendix C) on convection shows that the characteristic 
length for convection in nost flows 1S much greater than the 
conduction or diffusion lengths within the sheath and 
ambipolar regions. Consequently, convection may be entirely 
assigned to the ohmic study of electrode losses (see 
Appendix A). The boundary layer lies within the undisturbea 
plasma. That is, it lies beyond the ambipolar region, there 
is charge neutrality, and electron/fion concentrations are as 
would be predicted by the Saha equation{ 26] for plasmas at 


local equilibriun. 


In the modei presented here, the foilowing assumptions 
are made: 

1) Steady state, 
2) Chemical equilibrium in the boundary layer, but frozen 
flow in the sheath, 
3) No induced magnetic field, i.e., low Magnetic Reynolds 
number, 
4) Negligible ion slip, 


|) 3) a 
1 = 


6) No continuum radiation losses in the energy equation, 
7) No ion emiSsion, and 


8) Neutral and ion particle temperature is T , unaffected 
re) 


by Joule heating and uniform within the domain of the 
Sheath. 


B. CONTROLLING EQUATIONS 


1. Basic Equations 
For some MHD plasmas the mean-free-path 1s small 
enough that the sheath can be treated as a continuum. A set. 
of collisional equations is then used to describe the 


Sheath, ambipolar, and adjacent free stream regions. No 
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matching between these regions iS necessary Since they are 
generated ina self-consistent way within the computational 
eer ady . The equations used are a combination of the 
conservation equations and Maxwell's equations. 
Specifically, they are expressions of conservation or 


Charge, somentum and energy, as well as Gauss! law. 


A special representation of the electron nonmentun 
equation is embodied in a generaiized Ohm's law and may be 
Meeatten [| Ref. 26, p. 363] 


8 
— = — S i q 
Je Ugn.eVo TB] (J,xB) + D evn, (4) 


where the + sign applies to the electrons and - is for ions. 
The first term on the right-hand-side represents conduction, 
the second is due to the magnetic JXB current and the third 


term fmepresentsSs diffusion. The term Bb is the dHall 
S 


parameter and gives the relative effect of the magnetic 
Preld in the plasma. The eguation for species continuity 
can be written as 


ee (2) 
V re ne 


Aprendix cC shows that the characteristic length home 
1OnNization/recombination is much larger than the sheath or 
ambipolar region lengths. Under these conditions Fq. 2 


becomes 


Vie 0 {3) 


Poisson's eyuation comes from Gauss! law for the divergence 


Mian electric field and is given in potential form as 


wo = -(n,-n,) e/e 


) 


O 


Another eguation that will be useful later and which 
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applies under local equilibrium is the Einstein relaticn. 


Es 5 
Dee a/c (5) 


Electrons absorb energy at the rate of J *B per 
e 


unit volume from the electric field. Because of their 
relatively small mass, the electrons are inerficient in 
transfering energy to the other particles. For tnat reason 
much of this energy goes into raising the electron 
temperature above that of the heavier particles. Collisions 
With these particles, both elastic and ineiastic, tend to 
limit the energy retained by the electrons. fina quiescent 
plasma an éxpression accounting for the electron energy 
balance is (Ref. 26, p. 240-243) 





‘ee = S_e (6) 
Te E 3n k(T, T) L Vv 


2. Non-Dimensional Parameters 


The non-dimensional parameters developed in Appendix 


C are repeated here for use in this section. 


we ~ ni/N, nN, = n,/2, od = o/ >, 
= LY $= x/L f= y/L (7) 
where: 
6 = kT /e L = e-file ie and De = (beac Vee)? 
O O oO O ono 
T is the temperature of the neutral gas which is constant 
oO : 


throughout the region of interest, and independent of Joule 


heating. 
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Each of the three characteristic parameters oo», L 
O 


and non dpomerumettons Only Of a characteristic energy kT . 
O O 


fi now represents a non-dimensional electron charge, or a 
e 


measure of the deviation from the equilibrium density value 
in the free stream. f is a non-dimensional ion charge 


i. A i . e 
density, and within the sheath n #h , whereaS within the 


A A A A ne Cc. 
ambipolar region n =n. andno <n 7 and in the free stream 
e i e aha 
plasma fi = _— Mme aece sfOr purposes Of this work; the 
e i aha 
edge of the sheath will be where 1 approaches fi. within 


e at 

1%, and the edge of the ambipolar region will be where fh 

e 

approaches fi Wethin 2. 
Saha 


ee ee ae eee — 2 Se ee ee ee 


Using the non-dimensionalizing scheme on Poissuu's 
equation, Eg. 3 becomes 

A2 

V 


Sst? io. eee 8 
¢=n,-n, (8) 


=. ee ow oe mo ee 


Applying the non-dimensional parameters of Eq. 7 to 


E = -\/6 gives the following result: 


2 
igo Ge. an 
R= - —;—— Vd (2) 
e 


In the current equation, account must be taken of 


Srewettect of varying electron temperature. Since ver 
e iL 


the current is essentially an electron current. Neglecting 
magnetic field effects and using Einstein's relation, the 


Sirmrent if written as 
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rade 


= - ee = 
SS SS ea eT, Vn] (10) 
With the introduction of a dimensionless temperature, 


8 =T /T , the non~dimensionalized current equation hecomes 
e oOo u (kT ° 2 
e 

where 

1 me AA AA (12) 

7 nV¢ t evn, 
Equation 11 will contain the tern (1+/(32) in the aenoninator 
of the coefficient when magnetic field effects are included. 


miewaot product of Eqs. 9 and 11 yields 


yee S. 
Mercer). 6 ro 
e =} — Fo , ° (135) 
se wil J°V¢ 


henm—-dqinmensionalizing n and letting Coe 6 Me ave plOk, 
e S m 


w” 


mpomoantroducing Eq. 13 sas E@eese results in 
tL (ee ee 
e O (Oy ees 
aoe JY) = 
3ae 


a> 


_ (8-1) (14) 


As the electron temperature increases, the collision 
frequency with neutral particles increases affecting the 


terms [LL and @. Since electron-neutral collisions are by 
Qa 


: Mm 
far the most numerous (_ can be approximated by d= Vv 0 —@ 
eT 


ia 8kT 


a Oe ne?! 





Oo ey 2 
mm. ° ey 
e 
In the range of interest of the parameters, the ccllision 
cross-section changes little with temperature so @ depends 


He 2 
on 6 g 
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A modification to Compton's equation[{[27} for small 
values of E/P gives an electron mobility of the foro 
WZ 


ae es (16) 
€ eieton 
ee ITM, Cen™n 
Ve a: 
so that gces as 6 =) COtbinung Eqs. 4, 15, and 16 and 
ee 


using the perfect gas law results in the form of the energy 
equation used in this work. 
eo a la 

3°99 = -n, (87-8) a 


where A 9 

: 0.13(kT.) Eo mM, (18) 

7 26 2m 6 et 

P ven"eon 

Tiemeva Nee Or Sthe coerficient Y appears to dictate 

the degree of Joule heating in the system. This can best be 

seen by expanding Eq. 17 and expressing @ implicitely in the 

form 

9) Vo-Vo + 1 

= Ue (19) 

(y/6) Vo" Vn /n, + 1 

For physical reasons it 18S expected that @21. ic 

Same De yseen from Eq. 19 that for very small values of Y, 

6—>1, which is the constant temperature case. Se tage iv 


meee, Which is very large, an upper limit for 6 is found. 


—_ lla Pe a (20) 

Vo-9n, 
Miers interesting to note that for very large y, 98 becomes 
independent of Y, which means that it is independent of the 


chemistry of the plasma. 


Hewetnawe Y has been described, it will we useful to 


see what various MHD plasma compositions it represents. 
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Assuming atmospheric pressure, and the perfect gas law, iq. 
18 can be shown to be 
4 


A aaa 5 
-46 O i 1 


yY = 10 xX == 

Onen 

nen 
in MKS units where 7% is the atomic weight of the principal 
neutral species, or a maSS average of the neutral 


Sonstitucnts. 


For combustion type gases where the principal 


constituent might be es a typical value of Ycan be found. 


From Ref. 26, pages 139 and 148 for a temperature of 10900°K, 


fe) = 1.5x10-19me, 6) | on ei Seeniouid tou ne 
en n 


98.3. This type of gas has Y values which are low, 
characteristic of plasmas With high collision 


cross-sections. 


In many inert gas plasmas nitrogen is used to help 
the plasma reach equilibrium more guickly. The same 


reference gives typical values of Q = 5x10~20m2 and O = 
es n 


Muenen Y= 1.4x105. This is already at the upper limit 


described by Eq. 20. Even higher values are obtained with 
argon and cesium seed to the order of 10% because of the 
mameduer effect[ Ref. 27, p. 31]. 


— se ae ey fe eee Se aS eS SSS Se 


When Joule heating is small and collisions between 
electrons and neutral particles are sufficient for efficient 
energy exchange, electron temperatures will not differ 
pana ticantly Porn the neutrals. Furthermore, ion 
temperatures will be essentially the same as the neutrais 


regardless of Joule heating because of their large mass. 
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That is, the system of equations is stable even when T DT 


ESO 
Meeeause it has been assumed that T =T = a constant which is 
Ore 51 
unaffected Ey Joule heating. The following expression 
replaces the energy equation for this case: 
= = (27) 
TT, qT VB 


Now Eg. 1 1s solved for J andJ ina two-dimensional 





x iy 
CarteSian system with x and y coordinates. 
on 
= 2, ot a do S 
Jee = (itg”) UMS! Ox = 2) Ox 
a4 on 
= = 2 pS 
Sel ave tao “s= ayn (23) 
on 
Be Ai oo) S 
Jy = (1+8*~) [ Won.e me = De ay 
on 
oko S 
+ 8 (-u_new + Die )] 
S S's ox S ox (24) 


mieeretctucting Eqs. 23 and 24 into Eg. 3 and using Zqs. 5 and 
feeero nelp simplify results in 


on on 
*ng + Belge Be - = 


2 
eaey > - Vo-Vn, + (kT /e) V s ‘Fx By iy (By 


Equation 25 represents two equations, one for electrons 


(t+sign) and one for ions (-Sign). Now (3. = 0 because of the 
sf 


high 10n inertia. Hence the two eguations representing 
Spectres continuity and Ohm's law become 
on on 
, Cae 2 Pe a6 _ “Ne a6) _ 
nV ¢ Vo-Vvn, + (kT /e) V es B (xy ay sy ax) O (26) 


-n, V9 = Vo- Vn, = (kT /e) Von, = 0 (27) 


Non-dimensionalizing and substituting Eq. 8 into 26 and 27 


results in 


5). 





A A Pees On _ 
oe aet.) = 0¢ vn, + ) A, + BY = 0 (28) 
“w “Aw “Aw “Aw AN Don ms 
where Z 
on_ dag on_ dad 30 
OS see oe 
ox Oy, dy dx 
6. Species Equations with Non-Constant Electron 
Temperature 


For the case of significant Joule heating Egs. 28 
and 29 are nct descriptive of the species concentration 


Since T #T everywhere in the field. It can be shown that 


for this case the species equations become 





A Pe x AA Lats aan DESO) 1 A AA 
én, (n; n,) - oVn, Vo + 5 nV Vo + 6°V me 6VE hoe 
“aA A ~ on 
+e(- 5A, Se Ga, e+ oS (31) 
oy ox ox dy oy ox 
an an on 
ee os Me yd 938 2) 2 1 
ox OY ox oY oy dX 
A “A A AA AA ND) (32) 
n.(n.-n_) - Vn.-Vd - Von. = 0 
= 1 e cL 


Notice that for @ = 1 and\/9 = 0 these equations reduce to 
meeand 29. 


C. BOUNDARY CONDITIONS 


Many boundary conditions can be hypothesized for wall, 
electrode, free stream, and upstream/downstream Locations. 
Free stream conditions are actually those at the edge of the 
ambipolar region. They are taken to be zero potential and 


Zero space-charge, and a Saha-predicted total charge. The 
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electrode node potential 1s a fixed value. Other wail 
boundary conditions are not so easily determined since they 
are frequently vague and not well identified in the 
Meeebature. Also, since this model is new, the boundary 
conditions used must be compatible with the system that is 


being represented. 


Successful runs have been generated using the following 


@emartions with T =T: 
e O 


1) Free stream as mentioned above, 

2) Periodic electrode nodes allowing periodic conditions 
at upstream and downstream interfaces (Fig. 1), 

3) Line electrode nodes allowing current constrictions 
(the necessary condition for two dimensions [21}). Tne 
potential is fixed and ion and electron concentrations are 
zero at the node, as dictated by a catalytic surface, 

4) The inactive portion of the electrode wall is treated 
fmemeean “insulated wall". Current perpendicular to this 
wall is zero, and the wall eliminates space charge ‘ana 


stagnates charge motion. 


In reality a small sheath will form along the insulated 
Wall because of the difference between electron and ion 
mobilites. The analysis of Appendix C shors that the 
thickness of the sheath goes approximately as the square 
feet of the potential. Since the floating potential will 
generally be much less than the electrode node potential 
drop it is expected that the insulated wall sheath will be 
very much smalier in size than the electrode sheath. 
Conseguently, the hypothetical boundary conditions at the 


wall specified above are representative. 
A closer look at the electrode node boundary conditions 


is needed. Blue and Ingold[ 28] point out that though most 


authors use zero charge density at the electrode surface, 
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some cases require the modification n = 2J3J/({ec) for the 
s 


electrode poundary GOnaTtVens . Appendix G Give 


YU} 


moe 


characteristic values for J and © of 7x10% amps/m2 and 10 
c 


Wi 


m/sec, respectively, at typical temperatures. This resul 


(ae 
wt} 


in n= 10!9 particles/m3 which is of the same order as the 
Ss 


charge densities which will be encountered later in this 
WOrk. The difficulty in the application of this boundary 
condition comes from the fact that a true "point" electrode 
node would have an infinite current density. As will be 
seen later the numerical procedure uses a finite-~difference 
scheme and the electrode is not a true "point". The current 
density at that location becomes artificially dependent upon 
the mesh~size. To avoid this problem, either the charade 
density must be set to zero, or Sone new non-zero restraint 
must be found. The use of periodic boundary conditiois 
seems indicated since solutions of Lapiace's equation 
(Poisson's equation in the absence of a space charade) tend 
to be periodic in two dimensions. AS will be seen, several 
conditions were tried and the results are discussed in later 


sections. 


Other boundary conditions, such as a catalytic insviated 
wall were attempted with no reasonable results. This® is 


MmeeGussed £Lurther in the next section. 


Figure 2 is a representation of the domain in the 
Vicinity of the electrode, and some boundary conditions that 
are typical of those used in this’ work. The domain is 
represented numerically by a two-dimensional array of 
equally spaced points which are operated on by finite 


differencing the controlling equations. 


e e A A e 
Boundary conditions for o, De, sdna nh)  cenain the Same as 
e a 
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above for the Joule heating case. Temperature boundary 
conditions are straightforward. End conditions are again 
periodic and free stream requires that @ = 1. Since both 
the electrode node and walls consist of dense naterials with 
which electrons collide, it is assumed that the electrons 
give up their Surplus energy readily to these surfaces. 


Therefore, @ = 1 along that boundary. 
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Figure 2. Arrangement of electrodes and boundary conditions 
for computational array. 
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tie ae Poe OF SOLUTION 


— oy ey ee oe ee =e et eS ee ee eee 


moe LNTRODUCTION 


The governing eguations (namely Eqs. 8, 28, and 29) are 


e A a 
solved numerically for the parameters $, fi, and A. 
e i 


throughout the domain as Shown in Fig. 2. Later when the 
electron energy equation (Eq. 17) is included, a fourth 
parameter © is included in the field, with Eqs. 31 and 32 
representing the species equations. The method of sSoiution 
will now be discussed together with the application of the 
Meanadry conditions. The system of partial, non-linear, 
second-order differential eguations is baSically elliptic 
and constitutes a boundary value probilen. Tre partial 
differential equations are replaced by finite difference 
equations; the technigue for Solution is discussed telow. 
mie mon-linear nature of the controlling equations reyuires 
rather sophisticated computer techniques, and the boundary 


conditions present special problems. 


There are a number of possible approaches to the 
Numerical analysis of the equations and Severai were tried 
before a successful approach evolved. In order to better 
understand the nature of these complicated equations, this 
Section also discusses the unsuccessful appreaches and 


probable reasons for failure. 
Bee NON-LINEARITY CONSIDERATIONS 


The twce species equations contain the non-linear terms 
Which preSent special programming problems. One-step 
technigues for solving this system of three sinultaneous, 
secend order, non-linear, partial differential eqguztions are 


non-existent. One solution procedure in Such cases 18S to 
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include all non-linear terms on the "right hand side," that 
faye external to the coefficient matrix, and hope that these 
nonlinear terms change slowly enough with each iteration to 
render a convergent process. This procedure is known as the 
Jacobi method{ 29]. Looking at equations 28 and 29 for the 
@emetanit termperature case, or 31 and 32 for the Joule 


heating case, it can be Seen that only one term in each, 
es 2 2 7 e od 
either V n Sanya , 82s dinear, “and that there is no 
e i 


prescribed "load" on the right hand side. As will be _ seen 
in the expansion of these equations later, thLis leaves a 
dozen or more non-linear terms for the right hand side. 
Experience shows that the Jacobi method proves to be 


unstable. 


As a consequence of the failure of the Jacobi method, a 
quasi-Jacobi method 1s used. When the product £ two 
variables 1S encountered, one variable is treated as a 
constant coefficient for each iteration. This means that 
the non-linear terms are retained in the coefficient matrix. 
mae "constant" coefficients are updated alress every 
iteration, thus changing the coefficient fastest eee 
Successively higher voltages were calculated in. increments 


allowing smali changes to the system. 


ct 


~ 
J 


Had this guasi-Jacobi method not succeeded, the nex 
step would gave been to apply a Newton-Kaphsen technique 2 


3) 
to the equations. Because of the increase in the number of 
terms and the consequent increaSe in computer requirements, 


the Newton-Raphsen method was reserved as a last resort. 
C. GRID SIZE AND COMPUTATIONAL SEQUENCE 
The number of points reguired to define a domain 


sufficiently large to show the sheath and ambipolar regions 


1S critical (See the section on numerical results). hoe 
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example, if a two-dimensional field consisting of a square 
meeay Of 3% x 31 points is desired to be solved 
Simultaneously, then it can be shown that a finite 
difference scheme will require in excess of 66 million bytes 
Of storage for the coefficient matrix alone. Special 
storage techniques which take into account the  tkanded 
Structure of the coefficient matrix reduces the storage 
requirements to 12 million bytes. Even this is about 12 
times the core capacity of the IBM 360/67 at the Naval 
Postgraduate School. Obviously, to have a sufficiently 
large array and to keep within the core capacity the 
solution method must consider portions of the array at one 


tine. 


This problem is met by using the line iterative method, 
1.€@. by solving one line of the array at a tine, sweeping 
back and forth with successive field iterations. Diva @ den 
so, the Size of the coefficient matrix iS reduced to a small 
meet 2on of the full matrix. Because the egquaticns are 
elliptic, and therefore each point affects every other point 
in the field, it is desirable to reduce the computer time 
more 'anrormation™" to move through the field. This is 
accomplished by rotating the field 909° after every second 
field iteration such that the line sweeps go forward, then 
backward, then ud and then down. This approach takes 
considerably less execution time for a converged solution 
than the usual procedure without alternating directions and 


meration. 


The final grid size chosen for this line iterative 
method is 51 x 51. This decision was based upon the time 
required to reach a converged solution. The time for each 
run was of the order of four or five hours. The coefficient 
Matrix reguires only 51 x 3 x 10 x & = 12,240 bytes while 
the total space needed is 150,000 bytes for the constant 


temperature case, and 170,000 bytes for Joule heating. This 
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storage includes program space, variable storage and 


Pemutiron matrices for double precision arithmetic. 


To further increase the speed of convergence a method of 
successive over~relaxation is incorporated similar to that 
outlined in Ref. [30]. At the end of each line iteration 


the parameters are subjected to the equation 


= + ~ 33 
Acorr ee oad = 
where A is the corrected value of the parameter, A 2S 
Corr new 
the most recent calculation, A is the previous A Fars pees 
old COLE 


wis the over-relaxation paraneter. 
D. CONTROLLING EQUATIONS 


The finite difference representation of all terms was 
expressed to the order of h? when practical, where h is the 
grid spacing of the sguare 51 x 51 array. Ketter and 
Proweli (Ref. 29, p. 226-228) also describe the procedure 
mere setting up finite difference equations. Second 


derivatives in two-dimensions are given by 


2 
s = e e ry e ry e + e e - 4A, s h (34) 
VAS Bene) oil ,4 i, j4+1 eee i,j’ 


First derivatives have the forn 


oA 7 = 
x Pea 5 Pog a7 ae (35) 


f 


Where a first derivative is required at a boundary it is 


Meret en as either a forward or backward difference. 





an (36) 

= oe Vor » is ® aes iad je ° 2n 

Ox. . ( SE se De sg AE 

1,J 
Transforming Eq. 8 into finite agifferences gives 

na “A “A A “A Vin 

* ° ° ° e ° ¢ e ar e e + e 
Pi41,5 a *i-1,5 05 jt * 95 4-1 oh "hi ee 

- hn = 0 
e 


an 
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This is the simplest of the three controlling equations. 


The species equations are considerably more complicated. 
The finite dirference form of Eq. 28 is transformed here to 
illustrate the use of variables as "constant" coefficients 


fet ne non-linear terns. 


“An aN an an Nw nN Y 
(go. = ~ See {eA oul )-n ae at ean 
ty jt Ly Jot St A ese ore a aie eal 
w~ “aA “A Z. 
= poten ee a J78 
ay) FL TL yj) db Teas 
= on (n Sere (Otay ee Go: 4 a) In =n 
i i,j *i.3 ee he Hs 
* (\0 ea —n. )1/4n7 fa (n. +n. ) hn? (38) 
i,j+1 ie aye ee, 3 a 


This equation corresponds to a line in the y direction. The 
bracketed terms {} are treated as constant coefficients but 
are updated with each field iteration. The finite 
difference version of Eg. 29 resembles Eq. 38, but is 


considerably simpler since [3 #0. 
a 


When the electron temperature is uniform and equal to 
the neutral temperature, no energy equation per se is 


required since the condition T = T is incorporated in the 
eG O 


two species equations. However, for the non-eguilibriun 
case, the energy eguation is solved separately from the 
other tnaree. Hitters ses oO, ol, ahd s2 ane soived 
Simultaneously for assumed values of 9. With each iteratioa 


6 1s updated by solving Eq. 17. 


E. BOUNDARY CONDITIONS 
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Free stream boundary conditions are the Simplest. For 


A A A 
S@emmcas- the free Stream values are set at @= 0, nr =n = 
e a 
A 
n Peand © = 1. 
Saha 


Upstream/downstream conditions are chosen uO pe 
periodic. They require that not only the vaines of the 
unknowns pe the same at periodic points, but also tnat their 
derivatives be the same. For example, if the far upstrean 
station is labeled "1" and the far downstream station is "n" 


the two eyuations, then 


y 

tI 
ty 
o 
>) 
Cr 


(39) 
(40) 


memideapply for all parameters $, fh, fi, and 6, 
ee: 


The insulated wall is hypothesized to Le neutral and to 
Sieerate the charges, therefore, the conditions aire 


A fr n N ‘ E p ; 
Memem@ecadiiy set ton =n =A eC mW tl cod. t Lom ars 
e uh Saha 


taken from the restriction that perpendicular current is 


Oo 


Zero. From a non-dimensional form of Eg. 24 for zer: 


Surrent in the y-direction the boundary condition at tue 


wall is : 
_ o¢ an _ a¢ on | 
 —< +t — - Bi eter 2 irene: =O (41) 
oY oy Ox ox 
The x-derivative of f is zero along the wall andi = 
e e 
n So 
Saha es 7 - 
n 7 ps) fev Bhi an ay (42) 
A saha ~~ 
meena Chi ox 


The finite difference expression for Eq. 42 is 
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Nsaha?i,j 7 2 alpate spe 7 eeahe : ey 
year st (43) 


ay 


“NM sanali ,4+2 


“w 


+n + 
Seach 


BN aha ‘?it1, 5 Pi-1,4) 
The boundary conditions at the electrode node consist of 


Bree three conditions $ = } P n, = 0, and nm = C consistent 
fe) 1 e 


with catalytic electrodes. 


Pemeee OLUTION OF THE MODEL AND PRIOR ATTEMPTS 


Meemetactors affecting the choice of the calculation 
procedure and computational model are many. Boundary 
conditions, array size, electrode piacement, choice of 
coordinate system, and the form of the controlling equations 
meemeout a few of the things that had to be considered in 
modeling the problen. The final successful technique 
evolved by trial and error. The following paragraphs give a 
brief account of the evolution including successes and 


mar lures. 


At the outset a two-dimensional, flat plate, continuous 
electrode model WaS proposed. Because of the small 
thickness of the sheath relative to the electrode size, nost 
areas of the sheath wouid not experience end effects, and 
would be effectively one-dimensional. Further investigation 
led to the proof of the non-existence of a one-dimensional 


Solution presented in Appendix B. 


It then became obvious that the key to this problem was 
in two- or three-dimensional current constriction since the 
one-dimensional Cartesian geonetry offers no natural 
geometric means to decrease current density away from tne 


electrode. The electrode is then represented by a series of 


1) 


nodes aS shown in Fig. 1. An examination of the litecsaturc 
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showed that this mode of current constriction was Known and 
even helps explain the electrode spots which may be present 


on the anode. 


At first the model consisted of both the anode and 
cathode on opposing sides of a field vuf computational 
points, but it was realized that a field of this size couid 
mot possibly show all the details of the sheath. 
Conversely, a field capable of showing the required details 
could not contain both the electrodes. So the final model 
consisted of a series of line electrodes (for simplicity, 
the anode) along an “insulated” wall in a periodic field. 
Tae opposite boundary would then represent a free strean 
e@nmailton. The other two boundaries are upstrean/downstreain 


boundaries. 


Once the model was chosen, the job of solving the 
equations was undertaken. Initially, the three controlling 
equations (the energy equation had not yet heen considered) 


memeemeerved Separately for ¢, n , and n, respectively, in 


e a’ 
iterative fashion. The equations were Eqs. 3, 26 and 27 
With a Hall parameter of zero. The caiculations diverged 


almost immediately because in the vicinity of the electrode 
NMany derivatives are large and values of the parameters 


changed too rapidly from the initial "guessed" values. 


It became apparent then that the three equations would 
need to be solved simultaneously. Section IIi.C explains 
some of the problems of such a technique in terms of overail 
storage reguirements. i sero giteeore tits, Lio ekine 


iterative method was chosen. 
Armed With a Simple set of variables and the 


line-iterative method a solution was attempted. Convergence 


Was achieved, but very slowly, since it took many iterations 
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Pome iniLormation"™ to cross the array. So the array was 


rotated 90° at intervals causing the effective sweep to be 


back and forth, then up and down. 


At one time a coordinate transformation was effected 
putting three of the boundaries at eile de ea {The 
periodicity of the nodes was not yet considered.) This yvave 
no satisfactory solution, and the reason is believed to be 
the follcwing: PerevOoudias hd. 6 1S POoLSSon's equaticn, it 


Beauces to Laplace's equation in the ambipolar and free 


stream regions. Since there 1s no solution in 
two-dimensions for Laplace's equation Satisfying che 
requirement that ® = constant at infinity, then no solution 


to this problem was forthcoming. 


A variable mesh size was also attempted with some 
Success, but Since it was determined that the ambipolar 
region was of the same order of magnitude as the sheath, it 


WwaS rejected aS an unnecessary complication. 


Various boundary conditions were tried, but most cither 
did not work or gave physically meaningless results. Tue 
most difficult boundary to model is the insulated wall. 
From the final boundary conditions chosen the system of 


equations describes a certain phySical case. 


Once the present technigue was successful for various 
voltages, it proved its versatility by accepting the 
addition of a magnetic field and the energy eguation 
routinely. While not being the ultimate in possible 
technigues for soiving this system of equations, the present 
scheme has rendered useful results and promises to be 


further productive in other cases. 
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IV. RESULTS 


meee hOoCEDURE AND CONVERGENCE CRITERIA 


The simplest case, which included no magnetic field 
effects and negligible Joule heating began with the 
following assumed values: 

1. zero potential everywhere except at the electrcde, 
2. selected non-dimensional potential at the electrode, 
ae Zero charge density everywhere except at the 
boundaries, and 
4. Saha charge density at the boundaries. 
Using these assumed values, iterations for improved values 


proceeded. 


ju 


Periodicaily convergence was checked after many fielé 


iterations by comparing the left hand side of Eq. 8 with the 


} 


Mmegntehand Side at various-locations in the field. This 
convergence criterion may appear somewhat arbitrary, but was 
ereeenrwhen it was found that Eq. 8 was more @Gifficult to 
meewomy than Eqs. 28 and 29 and therefore represented 
Stricter convergence criterion. Additaoualliy, eis § 
contains all three unknowns providing a check of tail 


parameters. 


Mim error was calculated from the finite difference 


meeresentaticn of Eq. 37. 


E _ 9i,5+1 7 94, 5-2 


{I 
[—! 


h“(n ra) (44) 


As the calculations converged the error decreased with 
Subsequent iterations, but the rate of convergence slowed as 


the number of iterations became large. Eventuatly a point 


4g 





WwaS reached where further iterations were not cost effective 
as the additional improvement in the solution required 
excessive computer time. A second convergence criterion was 
used in addition to the the first. Certain potential values 


at critical points in the array were plotted witn successive 


field iterations. iets these values approached some 
assymptotic values, then the procedure was considered 
convergent. 

As discussed earlier, a successive over~relaxation 


technigue was used to increase the rate of convergence. It 
was found that over-relaxation parameters exceeding 1.71 
caused the field to diverge for all but very low voltages. 
Some improvement of convergence rate was noted for w = 1.05 
so this value was used for most runs. When attempting to 
achieve solutions for higher voltages or other conditions 
which might prove divergent, such as iarge Hall parameters 
or Joule heating parameters, a value of w= 0 was used in 
order to oput the least possible "stress" on the system of 


equations. 


Moen eGLIGIBLE JOULE HEATING AND NO MAGNETIC FIELD EFFECTS 


Convergence for this set of runs was slow. In order to 
achieve an error of less than 20% (f{€J<.2) over 1000 
iterations were reguired, even with a near optinum 
over~relaxation parameter of 1.05. This equates to about 
meur hours of computer on the IBM 360/67. in ea dd isetona) 
Meour (250 iterations) reduced the error to 16%. Final 
results shown in Figs. 3-18 were completed with [JE |[<.16. 
Further improvement would have required excessive computer 


time. 


meegees oa, 6, 9, and 12 are plots for the case } case 302 
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fies Valde represents one volt at T =2000°9K. Pictures 4, 7, 


re) 
Moree and 13 are for $ Seen Gn lgSao5m Go, Th, and 14 are 
meee = 29.02, mopGesenting 3 and 5 volts, respectively, at 
T = 2000°9K. 
O 


The overall size of the array for each of these runs is 
Mueeetimes the characteristic length defined in Eq. 7. At 
2000°K this neans that the dimensions or the sguvare field 
shown an the figures are 5.25x10-5m on a side. The 
non-dimensional charge density as determined by the boundary 
@e@mdi tion 1s a0-3. This corresponds to @.5x1017 
Meetrcies/m>, which 1S typical of an equilibrium charge 


density for seeded plasmas at this temperature. 


1. Potential Distributions 

Figures 3-5 are all one-dimensional profiles of 
potential and charge density distributions along a $line 
perpendiculac to the electrcde wall. The abscissa shows 
distance frcm the electrode in terms of the characteristic 
length. Fegures 6-8 are two-dimensional contour plots of 
the potentials for the three voltages. The electrode is 
located at the bottom center where the contours concentrate. 
Except for the absolute magnitudes, there is little change 
noted in the potential plots for the three voltages. AS 
Might be expected, the steepest potential slope occurs near 
the electrode where the greatest electric field exists. At 
the free boundary the slope of the fotential is nearly 
Bonstant indicating a constant electric field. A constant 
electric field is a property of two opposing infinite 
flat-plate electrodes. The plasma spreads the effects oz a 
POint or line electrode so that the free strean behaves as 


if the electrodes were infinite plates. 
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2. Charged Particle Distributions 
Figures 9-11 demonstrate considerable change for the 
electron densities with increasing voltage as do Figs. 12-14 


momeLoOn densities. 


An examination of the charge density plots in Figs. 
3, 4, and 5 reveal several interesting facts. First, since 
the edge cf the sheath is the line of points where the 
charge densities are equal, the figures put the sheath 
Pength at 35%, 90%, and 100% of the total length of the 
field. In the following table, the characteristic sheath 


lergth as calculated from £g. C.7 for T = 20009K is shown 
O 


next to the computer generated sheath length from the 
extimate of the convergence of the charge density lines in 


the figures: 


¢ d Cilaigectorlstic computer 
oO oO 

Sheath length generated 
-5 -5 

Sesh) l-0s vod tso6 00x10 m Pao Ft 0 if 
-5 -5 

7 7 ae 30 eVOb tse ls tox 10 m 4. 72x16 2 
-5 <5 

29.02 Oey Vokes loo oxat6 m bee OM mn 


These offer order-of-magnitude comparisons at best, but the 
predictions of the characteristic sheath length with Eg. C.7 
are useful in estimating the size of tne field for the 


ccmputer progran. 


Notice that the charge density slopes on the right 
Srerags. 4 and 5 are significantly greater than zero. Pus 
means that the field was not sufficiently large to encompass 
the entire ambipolar region. Since the boundary conditions 
Memeethe free stream are imposed aS conStantsS they are 
"forced" conditions and probably reSult in error for these 
cases. Gveomecesult of this is the erratic behavior of tne 


charge density close to the electrode. AS will he seen 
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Meer, this behavior has little effect on the current 


Seeetribution as long aS it occurs in a confined region. 


It is evident that the observation of 100% of the 
Pretd for the sheath length in Fig. 5 is probably low since 
emewecOluvergence orf the two charge lines occurs artificially 
at the boundary. To prove this, and the other observations 
made above, a run was made for the $ = 29.02 case doubling 
the field length to 1000 times the characteristic length, or 
foeoax 10-5 mm. Figure 15 shows the dramatic result. The 
erratic behavior of the charge density iines near the 
electrode is gone and the slope is essentially zero near the 
free boundary. The sheath length is measured at apbouc 
feeeexO-> mm, 2.3 times the characteristic predicted value. 
It is interesting to note that this is about the same ratio 
Mmetie reasured/predicted ratio of the $= 5.89 case. 
Figures 16-18 show the two-dimensional plots for this 


las 


@e= 29.02 case. 


Been ON-CONSTANT ELECTRON TEMPERATURE 


With the addition of the energy eyuation the effect of 
Moure heating on the system was studied. The value Y = 100 
was introduced first to Simulate combustion gases. The 
results are shown in Figs. 19-22 for the cases $ = 5.00 mnalc 
mee 02. In spite of the earlier discussion showing that the 
Size of the field is inadequate for the higher poteatial 
case, it is used here to avoid the conplication of changing 


the electrode spacing when comparing different potentials. 


The most cbvious change caused by the heating of tre 
electrons is the changed slope of the potential indicating 
less shielding for the Joule heating case. This is further 
shown in Figs. 23 and 24 which represent $ = 5.80 and 


f= 103. The charge density profiles are relatively 
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ater, this behavior has little effect on the current 


Meer tbution aS long as it occurs in a confined region. 


It is evident that the observation of 100% of the 
Meroe: Or the sheath length in Fig. 5 is probably low since 
mies CONVergence or the two charge lines occurs artificially 
meethe boundary. To prove this, and the other observations 
made above, a run was made for the $ = 29.02 case doubling 
the field length to 1000 times the characteristic length, or 
foeok 10-5 Mm. Figure 15 shows the dramatic result. The 
erratic behavior of the charge density iines near the 
electrode is gone and the slope is essentially zero near the 
free boundary. The sheath length 1s meaSured at anout 
Beegaxt0O-> im, 2.3 times the characteristic predicted value. 
It 1S interesting to note that this is about the same ratio 
Memes fcasured/predicted ratio of the $= 5.89 case. 
Figures 16-18 show the two-dimensional plots for this 


w 


®@ = 29.02 case. 


PemeeON-CONSTANT ELECTRON TEMPERATURE 


With the addition of the energy equation the effect of 
Joule heating on the system was studied. The value Y = 100 
waS introduced first to simulate combustion gases. The 
results are shown in Figs. 19-22 for the cases > pore ae ous! 
mo. O02. In spite of the earlier discussion showing that the 
size of the field is inadequate for the higher potential 
case, it is used here to avoid the complication of changing 


the electrode spacing when comparing different potentials. 


The most cbvious change caused by the heating cf tne 
electrons is the changed slope of the potential indicating 
Wess shielding for the Joule heating case. This is further 
Shown in Figs. 23 and 24 which represent Coes meena 


m= 103. The charge density profiles are relatively 
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unaffected by the addition of high temperature electrons, in 


fact there is no observable change in the sheath length. 


Figure 25 is the electron temperature profile (T /T vs. 
e o 


distance from electrode) using the same cut as the previous 
one-dimensional plots. This plot shows that the najority of 
the Joule heating occurs very close to the electrode where 
the largest electric fields exist. There is no  noticible 
electron temperature rise in the ambipolar and free strean 
regions. This graph is typical of the temperature profiles 
for the various cases and differs only in the height of the 
peak near the electrode. This peak appears to be a function 


of the potential as well as the Joule heating parameter. 


From the numerical standpoint, the addition of Joule 
heating cffers a benefit not easily foreseen. Convergence 
of the system with the addition of the energy equation was 
accelerated such that errors of less than 2% were achieved 
mietess than 250 iterations. The reason for this 
improvement is not yet clear but its advantages are obvious. 
The inclusion of Joule heating is a better approxinaticn for 
most generators, and provides the bonus of reducing the tine 


of numerical convergence. 


Pest rECT OF EQUILIBRIUM DENSITY 


An interesting case is that of increasing the 
equilibrium charge density to much higher yalues by Changing 
the boundary value of the electron and ion charge densities. 
Presented here are the results of setting i Sy ae soe 
boundary. At a temperature of 20009K this corresponds to a 
charge density of 102! particles/m3. Though not a comnon 
ase, it serves aS an indicator of some of the numerical 
properties of the method, and a good comparison for the 


Sheath length prediction. 
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Figure 25. Electron Tgmperature Profile, ¢=5.80, Y=100, 
=O, n=1lO-~, array size=500L, 51x51 grid. 
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A 


At fT = 2000°K and for $ = 5.80 and h=1 at the 
O 


Doundary, the predicted sheath length from Eq. C.7 is 
meek 10-* tn. At this same temperature the characteristic 
memoten as defined in Eg. 7 is L = 1.0x107-7 nm. Figures 2€ 
and 27 show the case where the field length is 51 or 
Mmoeznt0-° m. AS expected the sheath length is about 10% of 
the total field. 


Phe dotted portion of Fig. 27 iS drawn to agree with 
previous arguments on the electrode boundary conditions. 
Actually, some of these higher density runs were made using 
a different boundary condition, namely one that attempted to 
impose a Gonersoctent total cllrrent through the field. 
Though erroneous, it was found that this particular boundary 


memattion has little effect on the rest of the field, and 


eis 


h 


Ps 


ena t SenditLons were actually determined by 


() 


ct 


h 


electrode in several of the one-dimensional plots appear te 


1) 


non-conducting wali surface. The "kinks" seen close to 


be a resoluticn feature of the numerics vhich tends to 
override the effects of the charge density houndary 


mematcion at the electrode. 


E. feet, 6©OFC6U[U6U6THEh6UNNUUNBER6UGOF)6UELEMENTS.—©<CGIN THE ARRAY AND 
PeeCTRODE POSITIONING 


Earlier trials used a 41x41 array instead of the 51x51 
array used later. In general the smaller arrays did iot 
encompass the entire ambipolar region as can be seen in 
Pegs. 28-30 which are for nh = tad 6 = 11.6, and Peareay 
Size = 20L, 25L, and 35L, respectively, in 41x41 arrays. 
“They seem to lack sufficient room to satisfy the equations 


even after varying the mesh Size. 


Figures 28-30 also show the changes that occur when the 


mesh size is changed. The sheath length appears to be 
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@emorstent with the scaling, but the ambipolar slope does 
not appear wo Change, indicating the inadequacy of 
meSh-Ssizing to control the size of the field when the number 


oe elements 1S too low. 


For a qualitative look at a multiple electrode plot, 
mease 31-33 are provided. This scheme allows a current 
density comparison of this run with one of lower potential 
and the necessary symmetry in the absence of a magnetic 


field. Current profiles will be considered later. 


SeenON ETIC FIELD EFFECTS 


The introduction of a magnetic field induces the nost 
visible, although not unexpected, effects of ali the 
phenomena discussed so far. AS might be anticipated, the 
Symmetry of the field ais lost as shown in Figs. 34-43. 
Notice however, that while the potential shows Significant 
distortion, the charge densities do not, at least for the 


A 


case of 9 = 5.80. 


The angle of incidence of the potential contours with 
the insulated wall is approximately the complezent of the 
Hall angle (arctangent B) aS can be seen in Fig. 35 (459) 
and Fig 43 (279). This is easily understood by solving Eq. 


24 with J =0 and considering only conduction effects. The 


y 
mesuit is 
36 3¢ ea as 
ay / Ox ie 


This equation says that at points where J =0 and diffusion 


Hem small, the slope of the potential lines is -/3. In the 
Vicinity of the electrode where diffusion plays a more 


Significant fart the slope begins to deviate. 
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As With the runs in the absence of a magnetic field the 
eifects of the presence of a point electrode are dampened as 
the free boundary 1S approached. The magnetic field 
produces little or no distortion to symmetry outside the 
Sheath. 


freee CURRENT EROPILES 


For each of the previous runs a complete two-dimensional 
current distribution can be made using the potential and 
charge density distributions. The current lines make an 
angle with the potentiai contour lines approximately egual 
to the Hall angle and would theoretically be exact except 
for the effects of diffusion. Figures 44 and 45S show the 
current streamlines in the field for $¢ = 5.80 and f= 0 and 
1, respectively. These were sketched from numerical data 
from the "Current Progran" which gives the current in vector 
morm.e The streanlines entering the crieid from the right are 
from the active site that lies periodically to the right of 


the one in the field. 


To determine the total current flowing through the 
System it 1S necessary to look some distance frem the 
electrode. Theoretically, for this model the density of the 
current will vary from some given value at the free strean 
to infinity at the electrode, however, as explained for the 
case of charge density a "smearing" takes place due to the 


numerical approximation. 


lneeterns of the amount of information available, Figs. 
memeenrough 48 are the most important in this work. They are 
the current-voltage diagrams for the charge densities 
f = 10-3 and f = 1. The dimensionless current density was 
taken one or two stations from the free boundary (in order 


£0 avoid numerical errors at the boundary) and represents 


V7 
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DISTANCE ee ELECTRODE xL7! 


Figure 4h, Current Density Contour Plot, b= =5) 730 Rye OF B= =O). 
n=10-3, array size=500L, FIs, Seidl, 
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200 100 0 100 200 
DISTANCE FROM ELECTRODE «x L7! 


Figure 45. Current Density Contour Plot, §=5.80,Y=0, B=1, 
n=lO-, array size=500L, 51x51 grid. 


We, 


the macroscopic current that would be opserved in a 
two-dimensional model. The non-dimensionalization of the 
current waS consistent with Eqs. 7 and 11 and contain the 
term (14+) - The conductivity was introduced to replace the 
mobility, and oconStant terms were collected to give the 
sinplified coefficient shown in the figures. AN 
interpretation of these results as well as examples of their 


use is given in Section V as conclusions. 
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Ni CONCLUSIONS AND RECOMMENDATIONS 


eee ee —_——a_ ae = eS ee ce ee eee es 


Peel HYSICAL CONCLUSIONS 


A successful model of the Sheath evolved from the 
asSumptions of steady state, frozen flow, and uniforg 
temperature for neutral and ion particles. Qualitatively, 
the results of the analysis used in this work offer much 
insight into voltage drops attributable to electreades in 
contact with an MHD plasma. Quantitatively, these results 
are useful i1f it 1S remembered that the model 1s 
two-dimensional and does not give the additional degree of 
freedon for current expansion afforded by three dimensions. 
A summary of the more baSic concluSions which were drawn 
from the resultS is presented here. They wili be explained 
individually later in tnis section. 

1. The sheath can be self-generated from a consistent set 
of equations with or without the use of an energy 
equation. 

Die Tne electrode node boundary condition for charged 
particle density has little effect upon the field, while 
the insulated wall boundary condition has a profound 
effect. 

Bi Gireent cCONStrictions are necessary at the electrode 
to satisfy the system of equations. 

4, The resulting current-voltage diagram has a curvature 
@omsastent with theoretical predictions for space-charge 
in a one-dimensional flow of current and experinertal 
evidence. 

ae fMemeedrrent density, f£0F a given potential drop, 
varies inversely as the node spacing. 

6. The conductivity of the plasma within the sheath is 
critical in determining the sheath voltage drop, and 


mods to increase this conductivity will resuit in a 


84 


decrease of both sneath and boundary layer drops. 


1. Sheath Formation and Cnarge Density Profiles 

The results have shown that a sheath develops 
consistent with theory. The thickness of the sheath varies 
with the potential and is approximately that predicted from 
theory on Debye shielding. In each of the resuiting plots 
presented in the previous section the size of the ambipolar 
region is essentiaily of the order of magnitude of the 
Sheath length. No quantitative conclusions have been drawn 
as to the relative size of the ambipolar region other than 
to notice that higher voltages showed the tendency of the 
ambipolar region to increase in length relative to the 
Sheath. It is very possibie that at higher voltages the 
ambipolar region may grow to about the size of the houndary 
layer thickness. But since the potential drop in this 
Beqion 1S low compared to that in the sheath it is not as 
Significant in this analysis. The addition of Joule heating 
@amq a magnetic field did not appear to affect the sheath 
length. For the parameters used in this investigation, the 
difference between unlike charge densities (Space charge) 


was about 10% of the equilibrium density. 


Scme of the resulting one-dimensional profiles shown 
in the previous section exhibit erratic behavior in the 
first few stations past the electrode node. This is 
believed to pe caused by numerical convergence proplems in 
the procedure arising from steep derivatives near the node 
and the extreme non-linearity of the equations. THES 
appears to have littie or no effect on the rest of the field 
fPrce, as can be shown, the current distributions derived 
from these runs are consistent with those from smoother 
runs. Furthermore, charge densities other than zero were 


Mested aS a boundary condition at the electrode with no 
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change in the overall results. It seems then, that other 
conditions must be more dominant in the behavior cf the 


Potential field. 


One very Significant factor appears to be the choice 
for the insulated boundary condition. For example, a 


catalytic wall was attempted by setting nh =n =O, but no 
1 e 


reasonable resultS were obtained. The biggest problem with 
the catalytic condition appeared to be the non-compatability 
Memethie wall condition with the free stream boundary 
conditions, specifically the free stream slopes did not 
approach zero regardless of the mesh size or characteristic 
Pength. Piety Gradients occurring at a catalytic wall 
occur in a small region when compared with the sheath, and 
violate the assumption of a continuous fluid. Numerically 
this means that the charge densities must go fros zero to 
equilibrium values within one grid space. Physically the 
Momrcy Of the field should not "see" this wall anyway. 50 
the conditions chosen, namely that of an eguilibriun wall, 
meee consistent with the rest of the field and yield 


meaningful results. 


The current-voltage diagrams in Figs. 46 through 48 
Show the behavior of the system for various conditions. 
Figure 46 shows that the non-dimensional current density is 
inversely proportional to the electrode spacing since the 
current density for the $ = 29.02 case with the array 
size=1000L is approximately half that with the array 
miee=500. This is shown even more vividly in Fig. 47 where 
miremintroduction of multiple electrodes with doubie the step 
Bze (keeping the electrode spacing constant) resulted ina 
Samiale current line. The significance of this deopendancy is 


that the voltage drop for a given current density bececmes a 
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function or the electrode node density, and is thus a 
problem which has inputs from surface physics as well as MHD 
conditions. It must be remembered that for the square 
feos) BSC in this work the electrode node spacings are of 
the same crder as the sheath length. The minimum Spacing 
must be of this order to avoid one-dimensional effects that 
Meera cFesult from further "crowding" of the active sites. 
Spacing the nodes much greater than a Sheath Length makes 
memputaticn more difficult since a much Largez array would 
be needed to describe the effects. Thus, this work models a 
minimum Sfacing resulting in a near maximum cnrrent density 
for a qiven potential drop. On an actual electrode, the 


Spacihg wculd he random ang dependent cn empirical factors. 


mune f3 1S the current-voltage diagram for the 
Mei Case using a 51x51 array. It was noticed that for a 
given potential, the current density for tke B= 2 line was 
meproxitiately twice that of the B= POS oo wre Using this 
resuit a B = 2 line was sketched in Fig. He Deere ie 
converged solutions were available for this line at 


A 


pe 10-4. 


timesantroduction of the magnetic field moves the 
meres to the right as shown in Fig. 46, and appears to cause 
a smaller potential drop for a given current. But recall 
that the current is already multiplied by the term (143?) 
Sema the effective reduction in conductivity, and this 


partially offsets the reduction in the potential drop. 


The introduction or the energy equation shows that 
Joule heating increases the current for a given voltage, but 
Only to a slight degree. This is because the region that is 
affected by the Joule heating is small and very close to the 


electrode, well within the sheath. 


In tne absence of a magnetic field the 
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current-voltage diagram tends to be concave upward, a fact 
that can be predicted with the following arguments: Cobine 
(Ref. 27, p. 129) shows that the high pressur>2 svace-charge 
current should go as V?/y3 in one dimension. If y is taken 
to be the sheath Length in Eq. C.7 then it can be shown that 
V goes as Je. Figures 46 and 47 are nearly parabolic in 
appearance, and have curvatures consistent With this 


analysis. 


The zero voltage intercept produces an interesting 
problem since it has different signs for Figs. 46 ard 47. 
It would be expected trom Eq. 1 that for zero potential 
difference in the absence of a magnetic fieid, tne electron 
current density would pe positive because of the diffusion 
term and the positive charge gradient for the existing 
boundary conditions. However, recall that the validity of 
Mpemecatalytic electrode depends upon the fact that it is 
Charged. In the absence of a charged electrode, the node 
Should avpear the same as the insulated wall, and no density 
gradient should exist. Therefore, no special] siaqnificance 
can be placed on this intercept Without further 


investigation of the boundary conditions. 


Reference [31] reports on the beneficiai etfect that 
increasing the seed fraction has on lowering the electrode 
voltage drops. From the behavior of the electron collision 
cross-sections with temperature, the optimun 
conductivity[32] is reached at seed fractions higher than in 
the core because of the lower total tempetature at the 
eiectrodes. Heretofore, the only benefit would seem to be a 
higher bcundary layer conductivity, but as shown in Figs. 
46-48, the local conductivity aiso affects the sicath drop. 


Memee, a higher value of O is expected to decrease the 
O 


voltage loss associated with both the sheath and ambipolar 


regions. 
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An exampie would be appropriate here to show the 
orders of magnitude of the quantities being dealt with. 
Suppose it was desired to find the potential drop across the 
Sheath in a plasma where the temperature is 1800°K, the 
Current density is 0.2 amps-cm-2, free-stream conductivity 
1s 0.10 mhos/mn, and the charge density is known to be of the 
order of 19028 PabtLeres/ ns. Equation 7 shows the 
characteristic charge density to be 7x1029 particles/m3. 
Thus Fig. 46 would be the appropriate current-voltage 
diagram with f=10-3. The non-dimensional current density in 
the absence of a magnetic field is 1.50x10-5 which gives a 
dimensionless potential drop of 32.0. This corresponds to 


4.9 volts at the given temperature. 


Argyropoulos et al{ 33] cate the use of 4a 
sophistocated computer solution for the boundary layer drep, 
comparing the results with an experiment on the AVCO-APL 
Channel. Their resultsS show the anode drop to be 101 voits 
due to the boundary tlayer. They do not consider sheath 


effects. Using these data, nanely, T mie 2C00°K, 
We 


O = 3.95x10*% amps/ne, i= 1.02 and a chemistry censisting of 
toluene and oxygen with cesium seed, the conductivity at the 
eyeectrode iS 1.14 mhos/m. Figure 46 shows a non-dimensional 
@imrent density of 4.307x10-5. Using the [3 = 1 line this 
gives Ge=93526 for an actual drop of 6.13 volts at the given 


temperature. 


Another comparison is from the work of Sonju and 
Teno[ 34}. For this experiment they used the same chemical 
Penstituents as the previous one. The conductivity at the 
wall temperature of 18009K is 0.30 mhos/m and the otner 
parameters are J = 1.6x10“ amps/me and (one 1.8. Figure We 


myves j = 1./x10-*, 6 = 121.3, and a sheath voltage drop of 
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foes volts. 


Caution must be exercised when Making these 
comparisons since several factors have not been considered 
yet. First, tendencies towards non-~equilibriun would 
increase the conductivity near the electrode resulting ina 
smaller potential drop. Secondly, nothing has heen said 
about the node density of the experiments as compared to the 
assumptions made in the numerical solutions, namely that the 
distance between nodes is of the order of the sheath length. 
Shculd the actual node density differ from this it will 
mmpect, the current~voltage characteristics. Thirdly, the 
use of a two-dimensional model rather than a 
three-dimensional one nay tend to predict different voltage 


drops. The following example illustrates these limitations. 


An experiment by Kessler and Eustis[6@] was conducted 
using ethanol in oxygen with 1% KOH. Nitrogen was 
Mmeeroguced for cooling such that the Space ratio was 0.5. 
Pee nis run T = 1685°K, T = 27009K, J=.75x10%amnps/nm<-, 

wall core 
iS = 1.5, and the electrode conductivity is 0.56 mhos/n. 
Pigure 46 gives j = 3.697x10-%, jo 090s ane n de) =0gtes | 
mots. fies as a higher drop than the measured 45 volts 
mee Should account for both the sheath and boundary layer 


@emtributions. 


The primary difference between the above examples 
mepeats to be the conductivity at the electrodes. Since the 
Sailibrium conductivity has such a strong dependence on 
temperature, the Pec ToOdemeoUpDC Batu: cm sus cenit teat 
Parameter in controlling the sheath losses. A very accurate 
comparison can only be made using a plot generated from the 
correct charge density and on precise knowledge of anode 


ten peratures. 
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fiieetoseecommon practice to predict the losses 


according to the empirical formula 
V=A + Bé ) (46) 


Miere (0 is the boundary layer thickness and A and B are 
Memerants to be determined for a particular experiment. 
Such a scheme might be useful here if A is the sheath drup 
apa 6 2S the boundary layer drop. Using the comparisons 
above and from Appendix A for the boundary layer droz, the 
formula wculd look like 

18.8 + 11906 

6.1 + 11100 


where 6 is measured in meters and VY is measured in volts. 


Sonju and Teno V 


if 


Argyropoulos et al V 


fee UMERTCAL CONCLUSIONS 


The numerical technigue is discussed in detail in 
Section III. Generally it consisted of a finite difference 
scheme for the solution of a system of three (later four) 
coupled, non-linear, partial, second-order dirferential 
equations. The line iterative method of sclution was used. 
This method allowed a smaller storage requirement, and a 
choice of a 51x51 array balanced the conputer tine to 
converge with a mesh size sufficient to show the detail or a 
two-dimensional sheath. Convergence required four to six 
hours of computer time on the IBM 360/67 computer for each 
potential and/or Hall parameter chosen. The addition of 
Joule heating reduced the convergence time to one hour 
because of the damping effect of the temperature term near 


mye electrode. 
The use of successive over-relaxyation to reduce the 


convergence time met with only marginal success. Because of 


humerical instabilites, an over-relaxation parameter greater 
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moan 1.1 usually caused the process to diverge. 
Consequently, a value of 1.05 was used resulting in less 


than a 10% improvement in time. 


In its present form, the computational procedure 
outiined in this work is limited as to the size of the 
parameters describing the system. For example, divergence 
resulted for non-dimensional potentials exceeding 35, Hali 
parameters greater than 2, and Joule heating parameters 
greater than 10*. It can be observed from Eq. 28 thnat tne 
Hall parameter weignts the relative’ size of the axial 
derivatives compared to the cross-~-channel qerivatives, 
giving rise to numerical instabilities for excessively large 
values of . Tne maximum potential attainable seems to pe 
limited by the size of the array since a value of ¢ = 35 was 
@etarnedad for a 51x51 array and only 29 was achieved for a 


44x41 array. 


feeeeeCOMUMENDATIONS FOR FURTHER STUDY 


Since the principal limitations to this work arises 
because ci the large storage and time reguired Oe 
calculations, much more could be accomplished with improvea 
numerical technigues. A search by the author for a packaged 
subroutine capable of solving the computational matrix more 
efficiently than the one used here proved futiie. Improved 
Bechnigues are needed which will accomodate larger fields of 


computation. 


The validity of these results is limited by the fact 
that the current density decreases in two dimensions rathec 
than three. A three-dimensionai scheme, admittedly a major 
task to develop, could render more realistic information on 
the effects of the sheath and on other electrode phenomena. 


One possible scheme is to model a three-dinensional system 


S172 


By USING an axi-Symmetric cylindrical geometry for the point 


electrodes. 
Additional physical effects involving changes in 
boundary conditions are possible. FOr exampie, boundary 


conditions could change the anode to a cathode by changing 
the sign cf the charge and making the electrode an emitter 
Meee lectrons.{[Ref. 27, p. 106) The addition of convection 
into the systen of equations (should ve be deemed 
appropriate) 1S expected to accelerate convergence because 
it would ease the requirements of current density decreasing 


away from the electrode. 


Other effects such as neutral and ion heating «may be 
incorporated into the SGuetcmcO ives t) dtc Metue ria 
Meera oiiities known to exist for such a case. Conduction of 
heat in the solid walls, ion slip, and frozen flow may be 
mooed. It is anticipated that frozen flow woulda decrease the 
importance of the boundary layer drop and increase the 


importance of the sheath drop. 


With each aad ite nad effect the model rises in 
complexity and the numerical scheme becomes more unwieldy; 
therefore, only the principal effects should be reflected in 


the description of the sheath. 
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HNIQUE FOR DETERMING THE BOUNDAL 


4 
=—=-=— —- a <_< SS eee a ee ee eS —_— << = ee a ewe 


I 
E DROP IN MHD GENERATORS 


As a complement to the two-dimensional sheath prcecblen, 
and in order to assist in understanding the overail loss 
problem, a one-dimensional boundary layer analysis was 
investigated. The result is a method for determining the 
voltage drop across the thermal boundary layer which is 


Sinpler than most methods in existence today. 


The technique balances sophistocation with ease of 
Peemetcation. Because of its lack of rigor, many simplifying 
assumptions have been made which render the results 
maeceurate for certain cases. These assumptions are 
detailed in the following paragraphs, and experinentai 


comparisons are shown. 


1. Geometry 

Figure A1 is a schematic of the voltage drops across a 
typical channel. the geometry assumes opposing electrodes, 
but doeS not preclude the possibility of diagonal 
construction. This analysis looks at only the ohmic portion 
of the voltage profile and assumes that the sheath losses 
can and must be determined separately. The boundary layer 
may be fully developed and uses the wall conditions as one 


meethe boundaries for ohmic analysis. 
2. Boundary Layer Drop 


es EE we nee 


Reconstructing and modifying a formulation by Resaf{ 35], 


94 





/ 
/ 


AWUSONYd ADVLIOA 


x 2 


| eccmeamme | Vimee Lia 


—Wi-Why 


‘atTrzoad o8eqpToa (9 ‘*8oTeue 
TeoTAqoeTe (q ‘susTUeYydoW SSOT JO suofTdat BuyMoYs Yo Ieys (ev 
*1OZeADUSS GHW ou ur sesso, eseRpoa Jo suotzequsssiads1z sesazayy L’°V eansty 





Jako; AsDpunog 





9109 





4aAoj AsDpunog 


UsDoUS 





oe, 





the average resistivity (0-1!) across the channel can be 


written as 


_ D ae D-§ 
-l, _ 1 aE eel @: Cc 
(go -) = 5! CN B| J (ao aya J ore dy 
0 0 6 
= We (A. 1) 
+ f et ey 
D=5 


where O is the thernal boundary Layer thickness and D is the 
channel width between opposing electrodes. Assuaing 
constant conductivity in the free stream region, symmetric 
boundary layers at opposing walls, and non-dimensionalizing 


the integration variable, the equation becones 


1 
ae, . 1 He ee 5.10: (h. 2) 
(o ) Sg +2 f Ale) 2g 
0 


wh eed 


D 
S 


Rosa {Ref. 32, p. 75] assumed a plasma where the gas 
properties vary only across the channel, and averaged Ohm's 


law. The resulting equations for current are 


~8 U = 5 A.3 
Wy a (E., + U x B) pales ( ) 
= = '¢ . 
a. 0 Ey BU, (Ae) 
= WG 
where G= O (“a7 - B2 and all terms have been averaged 
memoss the channel (y-direction). The x-direction is the 
Barection of flow. Eliminating J , assuming (3 is a 
X 
Monstant across the channel, solve for E to get 
bi 
Me 2 5) 
E = UB + 4" “) = §6E As 
E, = Us ere (Cae SS) Sess 2he ( 
The actual output voltage is 
% ni 
= 2 
= = = eee: + BE D (A.6) 
oe = J Edy = ED = UBD 7 Bo ee. 
0 
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The aintercept of the voltage profile with the channel wall 
(see figure A1) is given by 
J 


_ Z 
BE, P= UP + GT (1 + BD + ELD 


oO 


(A. 7) 


femerag the difference and substituting for G the voltage 


loss is 
a, 2 “=i on 
Soloss = (U,-U) BD - (1+8°) [o lo ) - 1] x (A. 8) 


Memos loss is Ohmic and is attributed to the two boundary 
Mayet celrerature drops. For a turbulent boundary layer the 


first term will be small compared to the second 
JD ~— 


aS (1 BH} ic are 


bl 20 ) a ae) | (h9} 


Smeepcondaguctivity Ex 


oo en ee ee ee ee ee ew Se —— ye. = Se owes oe == ao 


O, 
An expression will now be derived for G to be used in 


Meese ?. «Hurwitz, Kilb, and Sutton{ 36] have shown that for 
cases where the current is entirely cross-channel (Faraday 
mode), the conductivity can indeed be treated as a scaiar 
Since for most such channels the magnetic field is 
approximately uniform between the electrodes. For the case 
fmeeacalar conductivity we can writef 26 ],[ 37] 

o = noe*/ [me .n, | ) XO) | (A.19) 


The subscript k denotes the various species present in the 


plasma, and x =n /n. For non-reacting plasmas the sum of 
k k on 
the products 78 will be relatively independent of 


temperature. We will assume this to be approximately true 
lao Ceacting plasmas. Taking the ratio of the conductivity 


at a point to that at free stream, Eq. A.10 leads to 


o7 





te 


oO Meeean 
Gin eco eon 
— (A.11) 
e“eo''n0 
Assuming a locally JWaxwelitan velocity distributvon, 
quasi-equilibrium (allowing the use of the Saha relation), 


mmomeconstant pressure across the boundary layer, Eq. A.1i1| 


becomes 
O 
=< o 97/4 Qo (1/0-1) {(A.12) 
s 
where (@ = ---- and 6 = T/T . 
2kT O 


4. Constant Specific Heat 

Miethais point a description of the turbulent boundary 
layer is needed. For the purpose of this analysis we will 
use the 7/7/th-power law and Reynold's) anatiogyf{ 38] tO 
describe the temperature field. An extensive and detailed 
study of the heat transfer problem related to tnese high 


temperatures was conducted by Brin[ 39], but in the interest 


of Simplicity the Prandtl Number is assumed to be one. FOr 
wemctant C , 
Pp 
eae oe - 
Weyer (7 TEES 
5 ae Lle- Os 
So ; 
Uy — ) 
a (=) = a as Suey, (A. 74) 
(1 - @) 


Mere combining Eqs. A.2 and A.12 and changing variables 
Meeording to Eq. A.14#, the result is 


w= il ae $774 4(1/9-1) 1 °F" )° 
(C7) == E +25 f 6 e ——_—— dé (2.15) 
c G (130) , 


6 
_% | 


Mmgerng £q. A.15 into A.9 for the average resistivity gives 


She 


= 2 
Vol = J (1 +8 )6/G. 
an 


= a Tie thG 
i fo rene se aol onde = 1 ey 
(1-6) Ae my 


peeeeenon—-Constant Specific Heat 


— =) ae ee ee ee ee ee ee ee ee a —_—s mee ently eee | ee rw Se ae ee 


For a plasma with a non-constant C the 1/7th-power law 


p 
is 
le) ania 
ae: 
= (~——-—) (A.17) 
6 Ag h 
A change cf variables would then give 
fee 3_(1+8°)6/o 
bl y e 
Mt ; -7/4 _a(1/8-1) 6 
—H—.. f 26 e (h-h )° de - 1 (A. 18) 
fh)’ 6 “i 
O W W 


Obviously an enthalpy table or equation is needed to solve 
fgee A.18. However, Mollier diagrams are obtainabie, and 
Once the entnalpy relation is known, the method shown in 


Simo dnalysis retains its simplicity. 


Peegumce A2 2S a graph of the conductivity from Sq. A.11 
for Toluene and pure oxygen with 2% cesium seed at a core 
temperature of 2600°K. Plotted also is the conductivity 
from a more sophistocated program provided by AVCO Everett 
which also assumes quaSi-equiiibrium, but does not make the 
other simplifications done above. The agreement at 2600°K 
results from the “forced" condition in the simplified 


program. 


An examination of Eqs. A.16 and A.18 shows that tne 
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effect cf changing the integration variable was to average 
MilemeresiStivity across the boundary layer in "temperature 
Space" instead of "geometry space". In each case this change 
generated a weighting PUTS E LOR which Weights the 


Contribution of the resistivity to the integral according to 


the temperature. Other weighting functions derived here 
are: 
7 (9-9,,) & 

mermeconstant C ; (1 fee = 

p (1 059% 

’ (hss) 

menevariable C ; ne Ome ee = 

Pp h p oy seal) ats 
Other descriptions are possible, and often desirabie. 
Useful descriptions for laminar boundary lavers are 


available in Cramer and Paif{ 40] and Kerrebrockf 47}. 


Figure A3 shows the effect that the weighting function 
has on resistivity. PRes PeCSrStivity . (Beer procal or 
conductivity) is plotted first from the AVCO program and the 
approximate vaiue from Eq. A.12. Next, the reacting plasna 
weighting function times the differential Ao is plotted. 
Lastly, the product of the first two graphs shows that the 
effect 1s to reduce the error considerably, by weighting 
more heavily those values near the core. In fact the 


average resistivity from Fig. A3 (area under the right hand 


eedon) 1S 0.0254 for the AVCO progran, and 0.0244 for the 


Simplified program. Obviously, the effect of the houndary 
layer shape is to smooth out the errors of the approximate 


resistivity calculations. 


From Eqs. A.16 and A.18 a non-dimensional boundary layer 


voltage parameter can be identified. 


2 (A.19) 
ane Ure  B ) 6] 


so that 
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=f ONO eRe) aI) oS: 
; (A. 20) 


Paemge (0) 1S the weighting function appropriate to the 
boundary tayer description. Figures A4ta and Aub are plots 


of ee vs. Wall temperature and core temperature for a 


constant C gas and CeSium seed. The fermer uses the 
- 


1/7th-power law while the latter uses a 1/8th-power law. As 
might be expected the 1/8th-power law predicts less of a 
voitage loss because the region of lower temperature plasna 
Bemeeceninner. Figures A5a and A5b are similar graphs fer the 
combustion products of Toluene in pure oxygen and 2% cesiun 


seed added. 


7. Expezinental Comparisons 

As with all experimental predictions the accuracy of the 
results is dependent upon the accuracy of the input data as 
weli as the abiiity of the model to repreSent the actual 
Situation. The TOL Ow dmg exanples Pistia te the 
convenience of this method, but at the same time show the 


importance of an adequate nodel. 


Sonjyu and Tenof 34) report the results of an experimental 
Mmemeeon the Viking I generator. Fig. A6é 1S a reproduction or 
the voltage profile for one run. Because oz burner 
innefficiencies they found the SOnaduct Pris yi of the 
Toluene/oxygen/cesium plasma to ke between 12 and 16 mhos/n, 
about half the 33 mhos/m that would have been predicted for 
complete combustion[{ 42]. Using {§ =1.8, ad (= 1.6x10' 
amps/m@, (= 12-16 mhos/a, 5 sim {haif the channel width 
for fullv developed flow), T =25009K, T =1800°, The 

core wail 


yoltage drop parameter is 0.282. This equates to a voltage 


loss per boundary layer of 119 volts for@=i6, and 159 volts 
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Figure A.4 Non-dimensional voltage drop $ versus wall temperature and 
free stream temperature for Cesium seed, constant ee 
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Figure A,5 Non-dimensicnai voltage drop OD yersus wall temperature and 
free stream temperature for Cesium seed, varying oe. 
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Figure A.6 Observed transverse voltage distribution. 
(After Sonju and Teno””) 
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MoeO,-12. Extending the voltage profile to the channel wall 


in Fig. A6 shows a on Oa DOWEn VOU SVOlcs. 


poether example comes from Argyropoulos et aif{33] who 

use a sophistocated computer technique to compare theory 

with experiments. They arrived at a boundary layer 

potential drop at the anode of 101 volts for an experiment 

@aectne AVCO-APL channel. Using the same data, namely 

T = 26009K, T = 20000K, 23+ 1.02, O60 om, a 
cor c wall 


chemistry of toluene, oxygen and cesium, and the 1/7th-power 


law, Fig. A.5 gives a non-dimensional potential drop of 


feed. This corresponds to a predicted drop of 110.5 volts. 


The next example shows the importance of selecting a 
Mmeeper model for the boundary Layer. Fig. A7 1s from data 
of Kessler and Eustis[6]}. The parameters for the experiment 
Wenres fF = 2.7/5 amps, electrode surface area = 3.68 cm, 


Z = 1685°K, T =27000K, G.= 10 mhos/n, §=1.5. The gas 
SL ere core 


was Ethanol in oxygen with 1% KOH added to increase 
conductivity. Nitrogen was introduced for cooling such that 


the 7° ratio was 0.5. Using the 1/7th-power law ae was 


meummda to he Q.290. Kessler predicted the boundary layer 
thickness to be 1.2 cm. The calculated voltage drop is then 
8.5 volts falling well below the 46 volts shown in Fig. A?. 
mimemaitference is too great to be attributed to sheath 


phenomena alcne. 


ite geometry of this last case is the cause of the 
discrepency. The channel depth is 3.1 cm while the width 
between opposing electrodes is 10 cm. This means that the 
flow does not behave as between two flat plates since the 
fudewalis have a strong effect on the electrode wall 


boundary layer. The model used in thiS comparison is 
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(After Kessler and Eustis ) 
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Seviously inappropriate. svelihibalewenestre | (GSvese qs | Sitt/ a oles P/5)) 
contains information on the problem of wall interference in 
Meieecil cular cross-sections and points out the increase in 
memettOn and boundary layer size. Teed EOss, 0 he hoe om 
approximation COGwEEroS procedure is to assune pat is Bay 
developed flow (as was suggested by Reseck[{ 43] for this same 
Channei). This results in a voltage drop prediction of 35.5 
volts. Ressler attributes 10 volts to the sheath drop for 
mene Cathode. When combined with the 35.5 volts predicted 


moewe, the entire drop is accounted for more fully. 


8. Conclusions 
Equation A.20 represents a useful means of determining 
thermal boundary layer voltage losses as lona as the 
foilowing conditions are met: 
1) A function of enthalpy vs. temperature is available, 
fee £OrL the case of £Eq. A.16, the gas has a fairly 


constant C in the range of temperatures of the channel, 
Y 


2) the core conductivity is known or can be estimated with 
some degree of accuracy, 

3) Hall fparameter, core temperature and wall temperature 
are known, and 


4) a suitable boundary layer description is available. 


The simplicity of this approach makes it desirable for 
integration with a computer program which predicts the 
performance of an MHD generator. Such prograns usually 
already contain free stream and wall temperatures, core 
Conductivity, and Hall parameters. Since these equations 
are one-dimensional they allow calcuiations at stations 
along the channel, eliminating the need to determine 
properties at points in a two-dimensional field. A sample 
program applying this method at one Station iS given in the 


section "Computer Programs." 
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APPENDIX B 


a Oy a a ee ee ee ee cee — 2 = _—= ee ees ee we Se ee ee ee ee ee ee ee 


we ee ww ee = 


feraerst attempt at solving a complicated set of 
equations is usually to reduce the problem to une dimension. 
For the equations which describe the sheath this procedure 
seemS reasonable due to the small thickness of the sheath 
relative to the electrode dimensions. Such an attenpt fora 
constant temperature plasma with no NGM WAcenkieyn “Ope 
recombination will fail because of tne nonexistence of such 
meescolution. Results published in Rer. 21 are summarized 
here along with Sable p pie skrersod ae conclusions which were 
instrumental in shaping the course of the ttro-dinensional 


Work. 


ie Order to simplify the analysis a change of variables 
is made. 
P = n, —n, c= n, +n, (Bal) 
The following manipulations are made on Eqs. 26 und 27 in 
moe text: 
feeiagnetic field effects are dropped, 
2) The eguations are considered in only one-dimension, 
3) The equations are added to each other and surtracted, 
4) The change of variables in Eq. B.1 are applied. 


Mirewresults, including Eq. 8 are 


2 
ao = fp (B. 2} 
dy? 
a aig, 4, BCE (B. 3) 
ape Be) = 0 

dy 
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oe (E 5a) + 2 = 0 (B.4) 
If Eqs. 6.3 and B.4 above are integrated and = = - do/dy 
is introduced then 
ie - 7 (B.S) 
- TE + ac = C (B. 6) 
- EE + i = .C, (B.7) 


[mames not difficult to 


momecucrrent 


snow that, 


Cy = “Cy = J/ ed.) 


where J 
y 


eGeeron ditfusion coefficient. 


C andcC are non-zero constants. 
1 2, 


the 
approacnes zero away from the anode, 
Ths 


Because ambipolar 


@erOo at the fplaSma proper. as 


|’ approaches 
& approaches 
E approaches 


Eventually a distance y = y 
fe) 


every yo y 
O 
|TE| < c,/2 


Hence 
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ms the cornventional e€urrent 


Obviously, 


region 


would be reached 


PM Ene anpseice Of a Net 


(B.6) 

flux and D so) Lie 
e 

since: J) -exisocec- 


is quaSi-neutral, is 
reaching a value of 
y increases: 
0 
ee 
Di oo 
SO, eat LOL 
(B.9) 





ag 

a 2 C,/2 Des caere te aa vo (B. 10) 

pegs oy) tt (C572) eee) (B. 11) 
which does not approach a constant as required. 


Furthermore, since 


= - Bere 
dy C, + SE ( ) 
P 

a 2 “CC, + [Ely,) + C)/2 x (y-y)1 EB ; (B.13) 


for y > y , which eventually grows at least as fast as 
Cc 
C 


E + (y-y,)* (B. 14) 


ana doeS not approach zero as required. 


mrmese ichViOUS inconsistencies arise because of the 


constant Bc The necessary reguirement for a solution is 


that tne current density decrease away from the anode giving 
mere to current constrictions at the surface. Under these 


Soma. tions, >. would decrease at an appropriate rate so that 


Tl’ approaches zero and ‘a approaches eo Mechanisms which 


allow a proper variation of oS are diffusion in two or three 


dimensions and ionization/recombination. The former negates 
the applicability of a one-dimensional CarteSian solution 


While the latter eliminates the frozen flow assumption. 





SOESEMON BIE), ( 1@: 


pee eee eee eee ee ee ee eee ee wre = we eee ee ce a 
SO ee oo oe ee ee 


Equations 1-6 in the text are cumbersone in that they 
contain many physical parameters which in their present form 
Mise be carried along as the equations are overated on. Tt 
iS useful to conduct a dimensional analysis or these 
Pavations in an attempt to collect the coefficients and 
Variables in such a way as to simplify the appearance of the 
equations and better understand their nature. Aliso, hefore 
introducing the problem into the computer it is imperative 
momexnow the sizes of the quantities that are being dealt 
With. For this latter reason, a fractional analysis shovwirg 
mre relative significance of certain physical effects 1s 


aiso given in this appendix. 


MeeeDimensional Analysi 


in 


The method outlined by Langhaar{44}]1is utilized here to 
determine the non-dimensional parameters relevant to tne 
Bieath equaticns. The significant parameters which appear 


Mebas. i and 3 are Jy , Ny Ye Cr E » pee , 3. and B, 
S O s 


where y represents any coordinate length. An expansion and 
Manipulation of Eg. 1 with Eqs. 3 and 5 will show that the 


| 


terms D and {{ can be replaced by kf. Since k and ?f 
S S O O 
e 


always appear together, they can be expressed as a singl 
tern. Additionally, the Hall parameter /f3 is already 
non-dimensional and will not be included in this analysis. 


Mrs leaves 6, n , y, ©, € , and kT. In the MKS systen 
Ss O O 


these parameters are combinations of the units of mass (4), 


length(L), charge (gq) and time (t). The matrix showing the 
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Meeckicients of these units for each parameter is: 


Yy n L kes 


a 
_ ase 


2 a3 1 0 = 2 | 


g -1 | 0 | 0 1 | 2 | 0 | 

te ae 0 I 0) | 0 Zi 2 

i ee eee in Die 
At first glance it might be expected that there are 
6 - 4 = 2 dimensionless parameters (number of parameters 
Minus the number of units necessary to describe them). 
However, closer examination shows that the t row is not 
linearly independent of the M row. In fact the rank of this 
matrix is three not four. THUS. Gthenpe “arc 16 — 3°53 
independent dimensionless parameters (number of parameters 
Minus the rank of the matrix) which describe the systen of 


Cgquations. This resuit would also have arisen if k and T 
O 


had been considered separately. 


Disregarding the t row, the homogeneous linear algebraic 
equations whose ccefficients are the numbers in the rows of 


maewaimensicral matrix are 


ay - an 1 a¢ = 0 
- 7 _ (Cr) 
2a, 385 + ax 32.5 + 2a 0 
“ay +- a) + 2a. = Q 


Serving for a, ome and 7 in terms of the other three gives 
uy 








2) | te 6a, = 2a, 
ag = ~3a, “ as (C.2) 
ag = “a, = 3a5 4F as 
Now setting the values aaa e and a a che besitos: 
a =1, a =O, anda =-1. 
4 
Continuing by assigning a value of one successively to 
a and a with the other two of that set equal to zero, the 
solution matrix is generated. 
d n i e e noe 
O 
os, a ee Se Geet | eee 
a : t i t 
al 7 | 0 i is 7 Q = 4 | 
n orf) -4 | Or. 6 8 | -~3 | 
ee ee 
A { 
y 0 OS” ae 1 | 1 
ia Bom Ge. Ses Se 
The resulting non-dimensional parameters are then 
Z Coulee 
rn ed o e 3 s 
‘i ii (= 0 ig y 
kT. KT SE6 en (Cc. 3) 
These terms are all independent since $, fi and ¥ appear in 


only one of e 


: 


o 


eae 


acne. In numerical form they Look like 
- ee ene b 
= P70 = cel C.0 ula ) 
_ _ =12 ar: (G4) 
= n/n, = 9.26x10 To n 
= y/L = 4 76x10° Ig Ny 
y > O 


The sheath size can be estimated by a technique called 


fractional analysis or approximation theoryf{[45]. Bo. 3 is 


Tae = mi 7 A) C7 e (C5) 


wom Characteristic" parameters 6 , ue onda, sae factored 
re) S fe) 


Out of the equation for each variable. 


Po 42 aon A 
<2 oe vas 7 
ae Vo E (n; ~n,) (C.6) 
S 
where 9 and n are the maximum value of the potential and 
O O : : 


charge concentration respectively, and ) is a sheath 
S 


mumGactersitic Length. Since in the sheath all the 
: ‘ Ee : 
Memeeres are of order unity, for the equation to be 
2 
| | neA ; | 
Seeeeeeerent the term ~-=--~- must also be of order unity. 
O 


setting this term equal to one and solving for mn , ‘ene 


~ 


nn? 


result is 





r ~* (C.7) 


This gives an ordec of magnitude for the charact¢ristic 
Seied@en length. For example, for a sheath drop of 10 voits 
in a plasma ionized to 1048 electrons/n3 the characteristic 


Sheath length is 2.4xi075n. 
3. Characteristic Curren 
Biom Fq. 10 the conduction current is 


ii Hen Vo 
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MoOEkKing in One-dimension and identifying a characteristic 


Surrent J , 4 fractional analysis gives the equation 
oO 


Ho EEN 


Mertrce that in order to be consistent the characteristic 


length in the potential derivative must again be A, the 
S 


sheath length. Solver rol rs a Whe nen = nts ome i et 
O O Q 


m@(volt-sec)—t, and @ = 10, gives J = 6.8x10% amps/m?. The 
O O 


meme Lor {jf comes from Ey. 16 for electrons colliding with 
e 


°. molecules at 1000°9K. 


4. Diffusion Length 

The potential is a monotonic positive function which is 
Suitable for approximation theory, however the charge 
density profiles have many inflection points in the sheath 
and are not weil suited. Nevertneless, it is possitle to 
consider profiles beyond the sheath where characteristic 


lengths can te detined. 


iMtmaitffusion is to play a significant role jn current 


pPreaduction then 


x on (e210) 
Uy Me we oy 
Again, using Peace eloma ts: analysis, the EOL Ova 
approximation 1S made: 
Cal 
J olga ( ) 


ee, Sh, 

AoveXt, 

Solving for I with the same parameters as above gives F 
d 


mers. 10-7. Since L< it can be concluded that any roie 
ad s 


played by the mechanism of diffusion should be visible 


Within the skeath. 





5. Tonization/Recombination Length 
From Hinnov and Hirschberg[ 46 ] the three-body 


recombination rate of electrons and ions is given by 


2 (C.12) 


where 
-33 ,kT,-4.5 er, 
n 


= . 1 ictal 
ee Te e (ony 


USing the technique outlined previously, 


i. on 


moe © = MPT a (ee (C. 14) 
i at e Oo e 
mete fF is the characteristic recombination time. tet the 
if 


term pee 0 (ki 7e) ~* eon? be of the order of unity and 


Sorve for 7. 


ie . 
a) —~35 kT, -4.5, 2-1 7 
tT, = (5.6x10 °°) () (n,”) eS 


At T=10009K and n = 710158m-3, the characteristic tine 
| O 


w= 


Mmecones 2.9x10-#sec. 


ie Oocrder to change the characteristic time to a 
characteristic length the mean thermal velocity is used 
Since this iS the transport mechanism for ionization ane 


Tecombination. Thus 





= eke (C. 16) 
Be = To = Ty Uh i 
e 
Using the same conditions as ahove L = 5.7x102 mn. Since 


r 


MeeeA it can be concluded that ionization and recombination 
s 
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do not have time to occur within the sheath for the 
Semattions used in this analysis. Since the results of this 
work show that the ambipolar region is of the same order as 
the sneath length, ionization/recombination is not expected 


to play a part in that recion either. 


Peeconvection Tine 

At atmospheric pressure and temperatures around 10009K 
sonic velocity is of the order of 103 n/sec. Pwis sas Sine 
mepene OL velocity of the flows of MHD generators. A rough 
characteristic time relevant to the sheath couid be 
estaplished by comparing the flow velocity with the size of 


the sheath. 


= (Ge) 
ue Ne he 


emma > = 10 volts and n = 1018 m-3, JT = 2.4x1U-8sec. 
O ) € 


PepGraracteristic sheath time can be defined iron 


=} /c ( 
ve \3/ 
for comparison with 7. Using representative values for i 
£ Ss 
ec, | = 1.2%10-19 sec. Thus, any effects occuring from 


S 


convective phenomena sutch as mixing and boundary layer 


effects do not have sufficient time to affect the sheath. 
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The principal program used in this work is the one that 
Caiculates the potential and charge density distributions 
for the steath and ambipolar regions. ThLS program is 


explained in Part 1 of this section and reproduced there. 


Auxiliary programs are those that calcuiate the current 
@istribution and produce the outputs used in this vyork. 


mite2 SHOWS the auxiliary programs. 


Peet Ss 1S a reproduction of the program usei in tne 
Calculaticn of the boundary layer in Appendix A&A. A sample 
enthalpy table is provided in data POushiae for the 


Toluene/oxygen/fcesiun mixture. 


1. Sheath Pregram 

The sheath program begins with an assumed solution and 
Operates on the field of parameters to replace the assumed 
solution with a more accurate one. Theoretically, the more 
iterations that are used, the more accurate the result is. 
Am example of an initial assumed solution is P{IT,J)=0, 
Meee? Sl (i,J) =BZ, and TH(I,J)=1.0 throughout the fiel4. 
(see the program comment cards for definitions) These three 
arrayS are operated on and a more accurate solution 15 


stored for later analysis or further iterations. 
The following pages are a reprint of this program, and 


an explanation of the parameters and subroutines is 


included. 
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Several auxiliary programs were used to process the 
meotres Of the sheath program. The first caiculates the 
non-dimensicnal current from the potential, charge density 
and temperature profiles. The second program creates a 
Graphical display of the potential and charge distributicns 
along a cut perpéndicuiar to the electrode and through it. 
The third prcgram uses the Same cut to regenerate the 
electron temperature ina one-dimenSional plot. The fourth 
program iS a modification of one by Biblarz{£47] which 
produces a two-dimensional display or contour map of tne 


potential and charge density fields. 


These pregrams are reproduced in the following pages and 


anclude explanations of their parameters and use. 
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3. Houndary Layer Proqran 
The program which solves equation A.18 from Appendix A 
MemeovOwn in the following pages. It calculates the solution. 


to the boundary layer voltage drop for a varying C. ac 


therefore reguires an enthalpy table for the calculation, 
and such a table is provided in data format following the 
program. This Canter e Ula table is LOL a 
Toluene/oxygen/cesium mixture at temperatures from 1500°K to 
264 0°RK. An explanation on the use of the progran is 


included. 
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